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A Discontinuous Galerkin Method for Parabolic Problems 
with Modified /ip-Finite Approximation Technique 

H. Kaneko, G. J. W. Hou, and K. S. Bey 
Abstract 

A recent paper [1] is generalized to a case where the spatial region is taken in R 3 . 

The region is assumed to be a thin body, such as a panel on the wing or fuselage of 
an aerospace vehicle. The traditional h— as well as hp — finite element methods are 
applied to the surface defined in the x — y variables, while, through the thickness, 
the technique of the p — element is employed. Time and spatial discretization scheme 
developed in [1], based upon an assumption of certain weak singularity of ||uf|| 2 , is 
used to derive an optimal a priori error estimate for the current method. 

Key words: Discontinuous Galerkin Method, Parabolic Equations, Modified hp- Finite 
Element Method. 


1 Introduction 

In this paper, the discontinuous Galerkin method is applied to the following standard 
model problem of parabolic type: 

Find u such that 

u t (x, t) — A u(x, t) = f(x , £), x £ D, t £ 

u(x, t) = 0, x € <9f2, t £ R + , (1-1) 

u(x,0) = u 0 (x), x £ fi, 

where is a closed and bounded set in with boundary R+ = (0, oo), A u = 
d 2 u/dx 2 + d 2 /dy 2 + d 2 u/dz 2 , u t = du/dt , and the functions / and u 0 are given data. 

*This author is supported by NASA- Grant NAG- 1-0 1092 
^ThLs author is supported by NASA- Grant NAG-1-2300 
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The discontinuous Galerkin method is a robust finite element method that can deliver 
high-order numerical approximation using unstructured grids. In this paper, region Q is 
assumed to be a thin body in i? 3 , such as a panel on the wing or fuselage of an aerospace 
vehicle. The traditional ft— as well as ftp— finite element approximations are used in the 
x — y variables, whereas, the p— finite element method developed, e.g., in [5], [15], is used 
in the 2 variable which describes the region through the thickness. The application of the 
p— finite element method through the thickness of thin structure, as compared to applying 
the ft — or ftp— finite element discretization to all coordinate directions, enables us to 
avoid structuring elements in R 3 that are too thin to satisfy the required quasi-uniformity 
condition (e.g. see [7]) that is necessary to deliver stable numerical approximation. We 
are coining the term ‘modified ftp ’-finite element method, as it differs from the traditional 
ftp— finite element method which uses ft— and p— finite elements on the same domain where 
the ft— finite element method provides a refinement of the region and the p finite element 
provides an enrichment. In Section 2, approximation power of the modified ftp— finite 
element method will be investigated. In Section 3, the discontinuous Galerkin method 
with the modified ftp-finite element approximation technique is established. Discontinuity 
is in time variable and time discretization is based upon the degree of singularity of ||ut|| 2 . 
The traditional ft— finite element method is employed in time. The convergence analysis 
given in [9] will be used. The reader is also reminded of recently published important 
paper [16] by Schotzau and Schwab in which various time discretization techniques are 
discussed. For instance, an exponential convergence rate in time of p— finite element 
method is obtained in [16] despite the presence of singularity in the transient phase of 
the solution. Time discretization used there is geometric. Schotzau and Schwab’s result 
extends the results in [3] and [4] in which no exponential convergence rate is reported. 
Also they discuss the ft— finite element technique in time using a class of radical mesh and 
obtain the algebraic convergence rate which is optimal. The radical mesh was chosen by 
analyzing the incompatibility between initial and boundary data. The present authors [1] 
established a similar time discretization technique for the discontinuous Galerkin finite 
element method, ft— version in time, which was based upon the singularity of ||ut|| 2 * Using 
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this analysis, it is shown in [ 1 ] that the optimal algebraic convergence rate in time of the 
discontinuous Galerkin method can be obtained under more dispersed, therefore more 
computationally stable radical mesh than the mesh used in [16]. 

2 Approximation Power of Modified hp Elements 

Let cj C R 2 and T C R be convex regions. For simplicity, it is assumed that T = [— |] 

where d = |r|. For simplicity, the thickness, d, is assumed constant over the domain. The 
Sobolev space of order k defined on u x T is denoted by H k (u xT) with the norm 

IMIifc,u>xr = ||D“u||2, 

0<|a|<fc 

where for each multi-integer a = (a 1( 02 , Q3). we have let |a| = ai + q 2 + a 3 and 

3 M 

D dx^dx^dx® 3 

We note that the Sobolev norm reduces to the usual L 2 norm when k = 0 . In this section, 
a best possible error estimate is derived for approximating an element in H k (uj x T) by 
the finite element function spaces. Let K^ v denote the master triangular element defined 
by 

= {(£,rj) € R 2 : 0 < r? < (1 + £)V3 -1 < £ < 0 or 

0<r/<(l-Ov / 3 0<£< 1}. 

Let S P (K^ V ) denote the space of polynomials of degree < p on K^ n , -i.e., 

S P (K = span{Crf:i,j = 0, 1, . . . ,p ;i + j < p}. 

First, the shape functions for the master element K ^ are formed. To accomplish this, 
the barycentric coordinates are introduced via 

Ai = (1 - € - i?/V3)/2, A 2 = (1 + £ - T7/V3)/2, A 3 = ( 2 . 1 ) 

A,'s form a partition of unity and A, is identically equal to one at a vertex of K^ v and 
vanishes on the opposite side of K^ n . The hierarchical shape functions on consists 
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of internal as well as external functions. The normalized antiderivatives of the Legendre 
polynomials are defined by 


M ) = /I Pl{t)dt ' (2-2) 

Now, the external shape functions consist of 3 nodal shape functions 

N i (Z,ri) = \ i , i = 1,2,3, (2.3) 

and 3(p — 1 ) side shape functions AfJ^(f, 77 ), i = 1, . . . ,p — 1 , j = 1 , 2, 3. The index j 
indicates one of three sides of Noting that ^(±1) = 0, 

&fo) = J(l ~V 2 )<Pi(v), i = 1,2,3,... (2.4) 

where is a polynomial of degree i — 1. For instance, <p\(rj) = — \/6, ip 2 (r/) = — \/l0i7 
and ^> 3 ( 77 ) = ^3(1 — Sr/ 2 ), etc. The side shape functions are constructed as follows: 

Nr^rf) = A 2 A 3 ^(A 3 - A 2 ) 

Af ] (£,7?) = AjA,v>i(Ai - A 3 ), i = l,...,p-l, (2.5) 

A^t?) =A 1 A 2 ^i(A 2 -A 1 ). 

From (2.4) and (2.5), there are 3+3(p— 1 ) shape functions. As dim(S p (K^ri)) — ( p+1 ^ p+2 \ 
the remaining ■ p ~ 1 ^ p ~ 2 ^ - basis elements are constructed in terms of internal shape func- 
tions. Clearly, nontrivial internal shape functions on exists only if p > 3. For p = 3, 
the bubble function on below serves as the internal function; 

&*<„&«?) = AiA 2 A 3 - ^=((1 - -|) 2 - £ 2 )- (2.6) 

Moreover, the collection I P (K^ V ) of higher-order internal shape functions can be con- 
structed from 

P(K^) = {b K(n v: v € ST*(Kto} = {b K( J ® 5 p - 3 (A c ,,), p > 3. 

Let T&, /i ^ 0 , be a triangulation of w. let x (L j , L 2 , L 3 ) and y L 3 ) be 

the element mappings of the standard triangle to the Zth triangular element K l E 

e.g., the linear mapping onto K l with vertices {(x*,i/*)}f =1 , 

Q l x {LuL 2 ,L z ) = ±x\L u C&LuL 1 M)=t.fa- 

t=i i=i 
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The space of all polynomials of degree < p on K l is denoted by S p (K l ) and its basis can 
be formed from the shape functions of S P (K^ V ) described above by transforming them 
under Q l x and Q l y . The finite element space 5 p,i (w, T h ) is now defined. For u>, p > 0 and 
k > 0, 

S p,k (u>, T h ) = {«6 H k (u): u\ K € S”(K), K € T h }. (2.7) 

Assume that a triangulation {T/,}, h > 0. of u consists of and that hi = 

diam(K l h ), for l = 1, . . . , M(h). 

In the 2 -variable for through the thickness approximation, the local variable r is defined 
in the reference element [—1,1] and T is mapped onto the reference element by Q z , i.e., 

T = Q 2 ([— 1,1]), z = Q z (t). 

Clearly, Q z is a linear function defined by 

2 = Qz{t) = ^(1 - r)(-^) + ^(1 + r)^, r € [-1, 1] 

The Jacobian of Q z is constant 

dz d 
dr 2 

In this paper, the basis functions of P p ([— 1,1]) are taken to be the one-dimensional 
hierarchical shape functions. See [15] for a complete discussion of the basis elements used 
in the p and hp- finite element methods. 

For example, in approximating an element in H l [— 1, 1], with l = 0, fa(r) = Pi-i(r), 
1 < z < p-h 1, where Pi- \ is the Legendre polynomial of degree i — 1, form the hierarchical 
basis functions. With l = 1, the external (fa and fa) and internal (fa, i > 3) shape 
functions are defined by 


Mt) = V> M t ) = 1 2 C 

= (^f 2 ) 172 III Pi-2{t)dt, 3 < i < p + 1. 


Note that fa's form an orthogonal family with respect to the energy inner product (•,*)£> 

(V'.iV’j \)e = J ^ ip'iiWjitfdt = j ^ Pi(t)Pj(t)dt = Stj. 
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Also note that the internal shape functions satisfy 

tpi{± 1) = 0, for 3 < i < p + 1. 

For the case l — 2 and p > 3, the four nodal shape functions and the remaining p — 3 
internal shape functions given by 

^i( T ) = K 1 _ r ) 2 ( 1 + T )> M t ) = j(l - r) 2 ( 2 + t) 

V’sCt) = -|(l + r) 2 (l -r), ip 4 {r) = |(1 + t) 2 (2 - r) (2.9) 

V’i('r) = (^r 5 ) 172 /_i(r - r?)Pi- 3 (r/)dp, i = 5,...,p+l. 

In this case, the internal shape functions satisfy 

dP \b ’ 

fy±l) = 0, for 5 < i < p+ 1 and j = 0, 1. (210) 

dr 3 

The nodal basis functions, ipi,i = 1, 2, 3, 4, in (2.9) also satisfy three of the four conditions 
in (2.10). For example, using the shape functions in (2.8), any element u G L 2 [-i, 1] can 
be approximated by u p G P p ([— 1, 1]), in the form 

1 — T 1 -j- 7" 

u p ( t ) = -^- U (-!) + “(!) + H a^T). (2.11) 

Z Z i~3 

For approximating the solutions of parabohc problems with the homogeneous Dirichlet 
boundary condition, the first two terms will be dropped, as u(— 1) — u(l) = 0. A sequence 
of triangulations {T h } h>0 is called the quasiuniform mesh if 

7 - all fe > 0, (212) 

with h = maxifer^ diam(K), and for some 7 > 0. P p (r) denotes the space of all polyno- 
mials of degree < p defined on F. The following is proved by Babuska, Szabo and Katz in 
[5] . See also [6] by Babuska and Suri on a related discussion. Here flo denotes a bounded 
polygonal domain in R 2 . 

Theorem 2.1 Let u G H k (Q o)- Then there exists a sequence z p G P p (Q 0 ), p — 1,2, , 
such that, for any 0 < / < k, 

i|u - 2plko 0 < Cp'^WuWkfio, 
where C is independent of u and p. 
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The parameters k and l are not necessarily integral. Their proof relies heavily on the 
approximation power of the trigonometric polynomials. 

With / = 0 in Theorem 2.1 and using the usual duality argument, the results in 
Theorem 2.1 axe further improved by Babuska and Suri in [6] (theorem 2.9), (see also a 
series of papers by Gui and Babuska [13]), to the ftp-finite element setting as follows: 

Theorem 2.2 Let Th be a quasiuniform partition n/fio* Then for k>l,uE H k ( f2 0 ), 

v€SP^T h ) “ ^11^°°) - ChV P~ k IMU*(n 0 ) 

where u = min(fc,p + 1). 

The corresponding error estimate in the || • j|//fc(Q 0 ) is also available in [6]. 

h- version in the x — y surface variables: First, the h-finite element approximation 
is considered in the x — y variables. Let z = s(r) = |r be the linear transformation of 
[—1, 1] onto T. Now consider the problem of approximating a function u G H k (ui x T) by 
a function from the tensor product space S[(w) <g> P p ( T), where 

SJi M = S r h (oj, Th) = {u6 L 2 ( w): u\K G S r (K), K € T h }. (2.13) 

For error analysis of h- version of the finite element method, the space S p ' k (u,Th) defined 
in (2.7) is not necessary, and the space S^(ui) of lower dimension can be used to attain 
the optimal convergence rate. Let P£:H 2 (uj) — + SJ^(u) denote the interpolation projection 
defined by 

(P[u)(x, y) = 'jr, u(x \ , yl)<p l i(x, y), for all (x, y) € K l and u € H k {w), (2.14) 

1=1 

where Th is a triangulation of u> with K l € T h and {(x(, y-)}, r = i is a set of nodes on K l 
with y?|(Xj, j/j) = Sij. Also, denote by Q p :H k (T) —> F P (T) a projection defined by 

p+i 

(Q p u)(z) — 530 ^( 2 ), for all z G T, (2-15) 

»=i 
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where ^(z) = ^i(s _1 (^)) for each i > 1 where ^ are defined, e.g., in (2.8) or (2.9). 
Recall that for k = 1, the constants a\ and a 2 are known in case of parabolic problems 
with Dirichlet condition, and it is assumed that a i > 3, in (2.15) are determined by 

I H z ) - ajVjiztfdz = min f \ u(z) - £ &j$;(z)| 2 dz. (2.16) 

J r j= 3 bjenjr J=3 

From approximation theory [17], 

||/-^ r ||L 2 (n) = 0(/ 1 r ), (2.17) 

Also Q p : L 2 (T) — > F P (T), from being the orthogonal projection in the sense described in 
(2.16) and from Theorem 2.1 that 

\\I-Q P \H(r) = 0(p- k ). (2.18) 


Let 


II Ph ® Qph = l|F h r (8) Qp||z. 2 (a,xr) = sup ||(F h r ® Qp)u\\ L 2 ^ x r)- 

\\u\\ 2 -l 


For u 6 I/ 2 (o; x T), 


r 

(■ Ph®I)u{x,y,z ) = y)^u(x[, 2 /[,z)^(x,y), 

i 1=1 


and 


p+i 


(Ph®Q P )u(x,y,z) = XICE^ 1 -!/)} 0 ^)^)' 


J=1 * t=l 


where a 7 depends upon u and obtained according to (2.16). First, approximation order 


under L 2 operator norm of ®Q P for P£ ® I is established. 


Lemma 2.3 For P£:H k (u;) -> S£(u;), 0 < r < k, and Q p :H k { r) -> P p (T) defined 
respectively in (2.14) an d (2.15), 

\\P r h ®I-P r h ®Q P h<Cp- k , 


where C is independent of p. 
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Proof: 


II Ph®I-Ph® Qp || 2 S sup ||u |, 2=1 ||(P h r ® I - P h r <g> Q>|| 2 
= sup N|2=1 {£i /*/ / r | ELi u(xj, yj, z)v>j(x, y)- 
E£J EL, d{x, y)a^j(z) \dzdxdy } 1 / 2 

< su P I |u|, 2 = 1 {E, /*< / r I EL, <P l i(x,y){u{x l i, rf, z) - ££±{ a j 'i> j (z))\ 2 dzdxdy} 1 / 2 

< su P||u|| 2 =i ELi{E, Sk‘ Jr Ip*( x > y)\ 2 \u{x[, y \ , z) - aj#j(z)| 2 dzdxch /} 1/2 

by Minkowski inequality 

< suPj i uii 2 = 1 {rmax 1 < i <r{V ; / A - ; | 'pj(x,y)| 2 !u(x[,y|,z) - Ej=! a^^z^dzdxdy } 1 / 2 

< C sup|| u | b=1 max,<i< r {/ r |it(x{, y\,z) - Ej=l a j 'H j (z)\ 2 dz} 1/2 

where M = £, J A , |<pj(x, y)| 2 dxdy 

< Cp~ fc , by Theorem 2.1. 

□ 

Similarly, the following lemma will be useful. 

Lemma 2.4 Let H k (w) — » with 0 < r < k and Q p : H k (T) — » P p (r). Then 

\\I®Q P ~P T h ®Q v h<Ch r 

where C is independent of r. 
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Proof: 


l|7 ® Q p - P r h ® <3 p || 2 = sup||„|| 2=1 1|(/ ® Q p - PH ® Q p )u|| 2 
— su P||u|| 2 =i{^^'eTh Ik‘ fr I Sjii a j{ x >y)'&j{ z )~ 

E^}(Ej=i a i( x l y))^j(z)\ 2 dzdxdy} 1/2 

where Ejii Oj(x, ^^(z) is the best L 2 (T) approximation of u(x,y, •) 

< sup|| u|b=1 {Eif'er h /*< /r I E£i {<*;(*> ») - EU y, ')</?>(*, y)}^j{z)\ 2 dzdxdy} x / 2 

< su P||u|| 2 =i EjiiiE/c'er,, //ri /r M*, 2/) - £[=i ylM*, y)l 2 l^j(*)| 2 dzd*dy} 1/2 

< C{Eifi € r h /it' 2/) - ELi a j( x l yl)<Pi( x , y)?\dxdy} l/2 

< C7* r , provided that a, G H r (cj), 

where the last inequality follows from a well known result of the approximation power of 
piecewise polynomials [17]. □ 

Using Lemmas 2.1 and 2.2, we obtain the following theorem which provides an error 
estimate for approximating an element in H k {ui x T) by elements from S^u/) 0 P p (F). 
The result will be used in the next section when the formulation of error estimate of 
the modified h — p discontinuous Galerkin finite element method for approximating the 
solution of the parabolic problem (1.1) is established. 

Theorem 2.5 Let u G H k (uj x T). Then there exists u* G S£(u>) 0 P p (r) such that for 
0 < r < fc, p > 0, 

||«-ti*IU 2 MD = 0(/ l r + p- fc ). 

Proof: Define u* = (T’f <g> Q p )u. Then, using Lemmas 2.3 and 2.4 

||ti-u*|| 2 = ||« - (Ph ® Q P )uh 

= ||n - (7 ® Q p )u + (7 ® Q p )u - (P£ ® Q p )u|| 2 
< ||« - (7 ® Q p )u\\ 2 + || (7 ® Q p )u - (/£ ® Q p )u|| 2 

= 0(h r +p~ k ). 
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□ 


hp- version in the x — y surface variables: Now we incorporate the hp- version of 
approximation technique in the x~~y coordinates. The goal is to approximate a function 
u G H k ( u x T) from the tensor product space S Plyk (u T/ l )0P P2 (r) for nonnegative integers 
Pi and p 2 * Analysis is similar to the one given in Lemmas 2.3 and 2.4 and therefore is not 
given. Using Theorem 2.2, it can be seen easily that 

Theorem 2.6 Let u e H k {u) x T). Then there exists u* € S pi,fc (u;, T h ) 0 P^T), 

\\u-u*\\ L2{ ^ T ) = 0{h''p- l k +p- 2 k ), 

where v = min(A;,pi + 1) and h — max^^j^ diam(K), with T ^ a triangulation of u. 

Remark: Let N(p) = ^-LKp+ 2 ) ^ Note that numbers of the degrees of freedom of 
S puk (uj,T h ) and Pp^T) are M(h)N(pi) and N(p 2 ) respectively. Since a single element 
through the thickness is used because of the specific structural consideration in this pa- 
per, we can not expect the total error to decrease by letting the diameter h —+ 0, -i.e., by 
letting the size of surface elements decrease to 0. The second error term would quickly 
dominates the overall performance of approximation in that case. In order for both of 
the error terms in Theorem 2.6 to decrease consistently, note that N(p) = 0(p?) and 
h = 0(M(h)~ l ). Thus the number of surface elements M(h) and the corresponding 
degree pi of polynomials should be selected so as to maintain 

M(h)~ k N(pi)~% ~ JV(pa)-4. (2.10) 

Equation (2.10) not only describes the consistent error estimates between the two terms 
but also indicates the consistent workloads between the surface and the through the 
thickness approximations. 
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3 Discontinuous Galerkin Method 


In this section, the discontinuous Galerkin (DG) method for problem (1.1) is developed. 
The discontinuity is introduced in time, which allows computation to march forward in 
time. This, when compared with the standard continuous Galerkin method, presents an 
enormous saving in size of computation. The DG finite element method for parabolic 
partial differential equations was studied in a series of papers by Erikson, Johnson and 
Larsson, [8, 9, 10, 11, 12]. In these papers, the convergence in time of h-finite element 
DG method is established for solutions which are smooth. More specifically, when the 
solutions are approximated by polynomials of degree r, then the algebraic error estimate 
of 0(At r+1 ) as At — * 0 is obtained. However, in many parabolic partial equations, 
solutions exhibit singularities at t = 0 due to the initial conditions. In a recent paper 
[1], the present authors established a graded time discretization scheme that captures the 
transient solution to optimal precision. The graded time mesh is selected by assuming that 
|M| 2 is weakly singular. A similar study of the graded time meshes is reported recently by 
Schotzau and Schwab [16]. They derive a set of graded time partition points by considering 
an incompatibility between initial and boundary conditions. It is demonstrated in [1] that 
the time discretization based upon ||u t || 2 provides more relaxed distribution of partition 
points. The paper of Schotzau and Schwab goes on to describe the p-finite element in 
time and obtain an exponential convergence in spite of a singular transient phase of the 
solution. We will not discuss the p-finite element in time in this paper. It will be taken 
up in [2] in which the complete p-finite element for parabolic problems is discussed. 

We begin by recalling several results from [2] that are pertinent to the present paper. 
The following conditions will be assumed. Recall from Section 2 that Q = u x T. Let 
(fr, T, S) denote a finite element discretization satisfying 

1. h is a positive function in C^Q) such that 

I V M x )l < Af, for all x G Q and for some M > 0. 

2. T = {K} is a set of triangular subdomain of u with each triangular element having 
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diameter such that 


c l^K — 



for all K € T, 


and associated with the function h through 


Ci hjc < h(x) < C 2 hn, for all x € K, K € T, 


where Ci > 0, c 2 > 0. 

3. 5 is the set of all continuous functions on ft which are polynomials of order r in 
x = (xj , x 2 ) on each K £ T and vanish on duj as well as which are polynomials of 
order p in the ^-variable in T. 


For the discontinuous Galerkin method for (1.1), we partition R + as 0 = to < t\ < 
■■■ < t n < ■ ■ ■ where we let /„ = with k n = t n — t., i . For each time interval, 

with q a nonnegative integer, We let 

K P = {v- R + - V hp : v\ In € P,(/„), n = 1, . . . JV}, 


where 



the space of all functions u* 6 S r h (ui) 0 P p ( r) 

or u* € S p ' k (u,T h ) 0 P p (r) such that > 

h = maxxeTh diam(K) where Th is a triangulation of ui 


and 

Pq(In) = {v(t) = Y, V i ti - V i € V hp}- 
i = 0 

The discontinuous Galerkin method is defined as follows: 


Find U such that for n = 1, 2, • • • , with ft = u 0 T, C/|q x/ „ € W% p and 

[ {(Ut,v) + a(U, v)}dt + ([[/]„_!, =f {f,v)dt for all v € W£ , (3.1) 

Jin J I n 

where [u>]„ = u>+ - w~, = lim s _ 0+ <-) w(t n + s), Uq = u Q , (u,v) = f n u(x)v(x)dx 

and a(u,v) = (VU, Vv). The smoothness of !|x 2 in) is subject to the initial condition 
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as well as to the boundary condition. For example, if we take in (1.1), uo(x) = n — x, 
f(x, t) = 0 and ft = (0, tt), then the actual solution of the corresponding problem is given 

by 

u(x,t) - 22 u °e~ j2t sin(jx), 
i= i 

where 

u j~n fo ( n ~ x ) sin (jx)dx 




= <>(})• 

In the following, C” s denote generic constants whose values change as they appear. Now, 

iMOiil = ikWiiLw = fce- 2 ^. (3.2) 

1=1 ai i= 1 

The last equality in (3.2) is justified because of the uniform convergence of YljLi Cj 2 e~ 2 ^ t 
with respect to t. Now using the fact that /” e~* 2 dx < oo, a simple change of variables 
(say, y = jy/2t) will show that the last expression in (3.2) is ~Ct~ 1/2 , which leads to 

IMOIfe = 

A similar argument shows that if Uj = O(p) for some initial value function u f) (x), then 
||ut(f )||2 = 0(t -1 / 4 ). This case arises when uq(x) = min(x, tt — x) for x € (0, n). If u° 
decays faster than j~ 2b as j — > oo, then l|u((i)[| 2 will be bounded as t — > 0. An initial 
phase for small t is the well known initial transient for parabolic problems. It is the case 
that the smoothness of the solutions of parabolic problems vary significantly in space 
and time with initial transients where highly oscillatory components of the solution are 
decaying rapidly. Therefore, in order for numerical methods for parabolic problems to 
be successful, it is imperative that the methods take a careful account of time and space 
discretization scheme so as to capture the transient solutions. An adaptive time step 
control scheme was established by Eriksson and Johnson in [9]. Time steps k n are defined 
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by controlling the size of 


min k 3 n \\uf ] \\i„ 

J<9+ 1 

where q is the order of spline used in time and = u t , = u tt , *4^ = and 

||t/;||j n = max t€ / n ||tc(^)|| 2 - Note that the method of Eriksson and Johnson requires some 
estimates concerning ||tx« and uf* = A u tt to achieve the second and the third order 
convergence in time. The approach given in [1] provides convergence of any order in time 
for the discontinuous Galerkin method by examining only the behavior of ||ut|| 2 - 

For 0 < a < 1 and q a nonnegative integer, define Q = For a positive integer N 
and T > 0, define 

C = (^) Q , n = 0,l,...,iV 

and 

t n = t* n T. (3.3) 

We let /„ = (t„_i,£ n ], n = l,2,...,N. Let k n denote the length of /„ so that 

K = I( ^)« - (^ )0 m n= 

Note that 


hence 


7~l 1 

k n < (J [— — T by the mean value theorem, 


k n <C 


N q+i 


(3.4) 


where C is a constant independent of n. The solution u(x, t) of (1.1) is then approximated 
in t over each I n by a polynomial of degree q. For example, with q = 1, let I' n w denote 
the linear interpolatory projection of u; € Hq in time onto Wm, viz, 


P n w(x,t) = - w(x,t n - 1 ) + - - w(x,t n ), for each t € 

kn k n 

Note that I ‘ n . considered as an operator defined on Hq is bounded with respect to the 
norm || ■ where 

IKi)lloo,/„ =max|k(i)IUoo(n)- 

Itin 
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Since ft is assumed to be of bounded domain, I l n is bounded with respect to || • ||/ n also. 
As was the case with the L 2 projection, equals the identity on polynomials of degree 
< 1. Expanding u(x , t) in Taylor series with respect to t at t n to the first or to the second 
order, we obtain, respectively, for each n = 1, 2, . . . , JV, 


\\u-r n u\\ In < f \\u t (s)\\ 2 ds. (3.5) 

J I n 

Lemma 3.1 Let 0 < a < 1, q a nonnegative integer and T > 0, we assume that t n , 
n = 1, . . . , N are defined by (3.3). Then 

nax / , 
n<N J ln 

where C n is a constant independent on N . 


max / s a ds < C n 


1 

N <*+ 1? 


Lemma 3.2 Let t n and k n be defined by (3.3). Then 

(1 + log ^--) 1/2 < yft, for each n = 0 , 1, . . . , N. 
kn 

Lemma 3.2 is used to guarantee the stability of the discontinuous Galerkin method. In 
the remainder of this paper, we illustrate the current 'modified ’ hp- finite element method 
by assuming the h-ve rsion in the surface x — y variables using the linear splines. Also 
we illustrate the cases for constant as well as linear degree in time approximation. Let 
{{xi,yi)}fi x is the set of nodal points which are the interior vertices of K in TV Let 
<Pj be the linear spline basis element defined by y t ) — 5^, for i, j = 1 , ... M. The 
superscript Z used in (2.2) will be dropped. For application of higher order spline basis, 
more nodal points are required over each K. The solution u of (1.1) is approximated by 

(t> 0) 

u{x,y,z,t) ~ZjL 1 {u(x j ,y j ,-% y t){l- ±) + u(x j ,y j ,%,t)(±-±) 

+ ££3 <4(t)ik{z)}<p j (x, i /). 

Note that u(xj, yj, 0), for j = 1, . . . , M are known from the initial condition. Also, for 
t > 0, the boundary values u(xj,yj,^p^,t) are given. As u(x, t) = 0, for x € dCl, t € R + 
in (1.1), (3.6) simplifies to 

M p+1 

u(x,y,z,t) ~ 5Z5Z a i(^i( 2 ) <Pj(x,y). (3.7) 

j= 1 i=3 
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At each time level £„, we approximate u(x, y, z, t n ) — u(x, t) by 

M p+l 

U n = U n {x) = U(x,t n ) = ^2^2a’i(t„)ip i {z)tp j (x,y) y n = 0, 1, . . . , AT. (3.7) 

j=l i=3 

To start the DG finite element method, we first require a] (to) and they are determined 
from u 0 (x). More specifically, for each j = 1, . . . , Al, since u 0 (xj, y jt z, t 0 ) ~ U°(xj,yj, z,t 0 ) 
E^a a J i{to)ipi(z), M(p — 1) many a- (£ 0 ) are found by solving 

/ ^^k{z)^i{z)dz ■ a J k {t 0 ) = / A u(xj,yj,z,0)ipi(z)dz, for i = 3, . . . ,p + 1. 

k=3 J ~~2 J —2 

Now, equation (3.1) can be formulated as follows: 

For n = 1,2, . . . , N, given f/ n_1,_ , find U = £/|/„ € P q (I n ) such that 

[ [(U u v) -t- a(U , n)]d£ + = f (/, n)df + , n”' 1 ^) (3.8) 

•'/n Jin 

for all t; € P q (I n ) where C/ 0 ’ - = ^o* 

For a special case, consider q = 0, -i.e., constant in time. As U n = f/ n ” = [/ n-1,+ in this 
case, (3.8) reduces to 

f v )* j (3.9) 

J In 

for all u G Fo(^n) and n = 1,2, ... AT. With (3.7), (3.9) becomes for each n = 1, 2, . . . 

M p+l M p+l 

j— 1 i=3 j=l i=3 

for each a = 3, . . . , p + 1 ; 0 = 1, . . . , M. 

For q = 1, we let U\ In = $„(£)+ ^ z £f i ’P„(x) where E£s af J (t n )^i(z) <£j(x, y) 

and ECa a? J (tn)ft(z) <Pj(x,y). As C/*" 1 * = and U n ~ l ' + = $„_i + 

(3.8) becomes 


//„[£(®n> V )+ a($» + L ^= i ^„,t?)]d£ + ($„,«" J ) 

= //„(/) v )dt + ($n-l + ^n-ljf” -1 ) 


(3.10) 
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for all v € Pi(/„). By taking v — tp a fg and (3.10) reduces to the following 

linear equations for 2 M(p — 1) unknowns of J (f„) and af J (t„): 

Ejii Ef=3 af' J (t n ) { v>^«) + k n a((pjtjji, <fpil> a )} 

+ Ejil Ei2$ a f’ 3 ( t n){(<Pjl/>U <^a) + k n a(lfij 1 pi, Vgipa)} 

= //„(/(*)> V^a) + Ejil Ei^k^C^n-l) + WM> 

a = 3, . . . ,p + 1; /? — 1 , . . . ,M 

Ejil Ei =3 0 ,f’ 3 (tn) ! f a {(Pji’h <P 0 ll>a) + k n a(<PjXpi, <pf} 1 p a )} 

+ E^lEi’=3 a t' 3 (*n){ \ {<Pjll’i,<Pp1p a ) + %a(<P j 'l’ i ,<P0*l’a)} 

= K //„(* - *n-i V 0 lp a )dt, 

a = 3,...,p+ 1; 0 = 

The foiiowing theorem can be proved by minor modifications to the proof of theorem 
1.1, [9] and by making use of Theorem 2.5. The present theorem is described for fi = 

uj®t c r 3 . 

Theorem 3.3 Suppose that there is a constant 7 such that the time steps k n satisfy 
k n < 7 fc n+ i, n = 1 , . . . , N— 1 and let U n denote the solution of (3.8) approximating u att n . 
Here u is approximated by a polynomial of degree q > 0 over each I n for n = 1 , . . . , N — 1 , 
and u(*, •, *, t) is approximated by an element from S£(u;) 0 P p (T) for each t G R+, where 
u is a polygonal domain in R 2 . Then there is a constant C depending only on 7 and a 
constant ( 3 , where px > 0 hx and px is the diameter of the circle inscribed in K for all 
K G T h , such that for n — 1, 2, . . . , AT, 

\\u(t n )-U n \\ 2 < <7(l+log^) 1/2 {max||u-7>|| /m +/i 2 || J C>2 y u|| /n +p- fc ||u|| /ni//fc(r )}, (3.11) 

where |M|/ n ,tf*( r ) = max fG / n H^WIIn^r) an d D 2 y denotes the second order derivative with 
respect to x and y variables. 


Lemma 3.2 guarantees that the current DG finite element method with the graded tem- 
poral meshes defined in (3.3) is a stable scheme. Also Lemma 3.1 provides a bound for 
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the term rnax n < m \\u - I l n u\\ Im in (3.11) provided that ||u t || 2 = 0(t~ Q ) as t -> 0. Theorem 

2.5 is used to control the last two error terms in Theorem 3.3. In summary, we obtain the 
following theorem which utilizes the traditional h-finite element in the surface variables. 

Theorem 3.4 Suppose that u £ H k (u x T) and ||u t || 2 = 0(t~ a ) for 0 < a < 1 and that 
the time partition points t n are taken according to (3.3). Let U n denote the solution of 
(3.8) approximating u at t n . Also assume that u is approximated in time variable by a 
polynomial of degree q> 0 over each I n for n = 1, . . . , N—l, and w(-, •, £) is approximated 
by an element from Sl(u) ®P p (T) for each t £ i? + , where u) is a polygonal domain in R 2 . 
Then 

||«(«n) - U n \\ 2 = 0(N-W + h 2 + p~ k ). 

If higher order r > 2 splines are used in then the second term in the error can be 

replaced by h r provided ||D^ y u||/ n is bounded. 

In the case of the /ip-finite element approximation for the surface variables, Theorem 

2.6 is now used to establish the following. 

Theorem 3.5 Suppose that u € H k (uj x T) and ||ut|| 2 = 0(t ~ a ) for 0 < a < 1 and that 
the time partition points t n are taken according to (3.3). Let U n denote the solution of 
(3.8) approximating u at t n . Also assume that u is approximated in time variable by a 
polynomial of degree q > 0 over each /„ for n = 1, . . . , N—l, and u(-, -,-,t) is approximated 
by an element from S Pl ' k (uj, 7),) <8> Pp2 (T) for each t £ where u is a polygonal domain 
in R 2 . Then 

ll«(«») - U n \\ 2 = 0(N~ ^ + Vp^ + P f k ), 
where v — min(fc,pi + 1). 

Numerical experiments of the presently proposed ‘modified” h — p finite element 
method for parabolic equations will be reported elsewhere in future. 
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